3.2 Exploratory Analysis of clusters
Load R packages and define colour functions
library(tidyverse) ; library(reshape2) ; library(glue) ; library(plotly) ; library(dendextend)
library(RColorBrewer) ; library(viridis) ; require(gridExtra) ; library(colorspace)
library(GGally) ; library(ggpubr) ; library(ggExtra)
library(WGCNA)
library(expss)
library(polycor)
library(rstatix)
library(biomaRt) ; library(clusterProfiler) # For the ORA
library(foreach) ; library(doParallel)
SFARI_colour_hue = function(r) {
pal = c('#FF7631','#FFB100','#E8E328','#8CC83F','#62CCA6','#59B9C9','#b3b3b3','#808080','gray','#d9d9d9')[r]
}
# Get colors from the ggplot palette
gg_colour_hue = function(n) {
hues = seq(15, 375, length = n+1)
pal = hcl(h = hues, l = 65, c = 100)[1:n]
}
# Assign an HCL rainbow colour to each module
get_mod_colours = function(mods){
n = length(unique(mods))-1
set.seed(123) ; rand_order = sample(1:n)
mod_colors = c('white', gg_colour_hue(n)[rand_order])
names(mod_colors) = mods %>% table %>% names
return(mod_colors)
}Load preprocessed dataset (code in 2.1.Preprocessing_pipeline.Rmd)
# SFARI Genes
SFARI_genes = read_csv('./../../SFARI/Data/SFARI_genes_01-03-2020_w_ensembl_IDs.csv')
# Load Gandal dataset
load('./../Data/preprocessedData/preprocessed_data.RData')
datExpr = datExpr %>% data.frame
# Clusterings
clusterings = read_csv('./../Data/preprocessedData/clusters.csv')
# Updates genes_info with SFARI information and clusters
genes_info = genes_info %>% left_join(SFARI_genes, by = 'ID') %>%
mutate(gene.score = ifelse(is.na(`gene-score`) & Neuronal==0, 'Others',
ifelse(is.na(`gene-score`), 'Neuronal', `gene-score`))) %>%
mutate(Group = factor(ifelse(gene.score %in% c('Neuronal','Others'), gene.score, 'SFARI'),
levels = c('SFARI', 'Neuronal', 'Others'))) %>%
left_join(clusterings %>% dplyr::rename(Module = Modules_dh), by='ID') %>%
dplyr::select(ID, baseMean, log2FoldChange, shrunken_log2FoldChange, significant, Neuronal,
gene.score, Group, Module) %>%
mutate(module_number = Module %>% as.factor %>% as.numeric)
rm(clusterings, dds)3.2.1 Cluster-level metrics
Associations between clusters and sample traits
This relation is calculated using the correlation between each cluster’s eigengene and a vector with value for each sample of the trait that wants to be studied.
In the WGCNA documentation they use Pearson correlation, I think all of their variables were continuous. Since I have categorical variables I’m going to use the hetcor function, that calculates Pearson, polyserial or polychoric correlations depending on the type of variables involved.
I’m not sure how the
corPvalueStudentfunction calculates the p-values and I cannot find any documentation…Compared correlations using Pearson correlation and with hetcor and they are very similar, but a bit more extreme with hetcor. The same thing happens with the p-values.
datTraits = datMeta %>% dplyr::select(Diagnosis, Region, Sex, Age, PMI, RNAExtractionBatch) %>%
dplyr::rename('Batch' = RNAExtractionBatch)
# Recalculate MEs with color labels
ME_object = datExpr %>% t %>% moduleEigengenes(colors = genes_info$Module)
MEs = orderMEs(ME_object$eigengenes)
# Calculate correlation between eigengenes and the traits and their p-values
moduleTraitCor = MEs %>% apply(2, function(x) hetcor(x, datTraits)$correlations[1,-1]) %>% t
rownames(moduleTraitCor) = colnames(MEs)
colnames(moduleTraitCor) = colnames(datTraits)
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nrow(datExpr))
# Create text matrix for the Heatmap
textMatrix = paste0(signif(moduleTraitCor, 2), ' (', signif(moduleTraitPvalue, 1), ')')
dim(textMatrix) = dim(moduleTraitCor)
# In case there are any NAs
if(sum(!complete.cases(moduleTraitCor))>0){
print(paste0(sum(is.na(moduleTraitCor)),' correlation(s) could not be calculated'))
}
rm(ME_object)Modules have very strong correlations with Diagnosis with really small p-values and not much relation with anything else. Perhaps a little with PMI and Brain Region.
# Sort moduleTraitCor by Diagnosis
moduleTraitCor = moduleTraitCor[order(moduleTraitCor[,1], decreasing=TRUE),]
moduleTraitPvalue = moduleTraitPvalue[order(moduleTraitCor[,1], decreasing=TRUE),]
# Create text matrix for the Heatmap
textMatrix = paste0(signif(moduleTraitCor, 2), ' (', signif(moduleTraitPvalue, 1), ')')
dim(textMatrix) = dim(moduleTraitCor)
# There's one entry of the textMatrix that is not being rounded off properly, so I'll have to set it manually
# used to be 0.6 (8.99999999999999e-304)
textMatrix[14,1] = '0.6 (9e-304)'
# Change name of gray cluster to reflect it's not a cluster
yLabels = paste0('Cluster ', rownames(moduleTraitCor) %>% as.factor %>% as.numeric)
yLabels[yLabels == paste0('Cluster ', ncol(MEs))] = 'No cluster'
labeledHeatmap(Matrix = moduleTraitCor, xLabels = names(datTraits), yLabels = yLabels,
yColorWidth = 0, colors = diverging_hcl(21,'Tropic') %>% rev,
bg.lab.y = gsub('ME', '', rownames(moduleTraitCor)), xLabelsPosition = 'top', xLabelsAngle = 0,
textMatrix = textMatrix, setStdMargins = FALSE, cex.text = 0.8, cex.lab.y = 0.75, zlim = c(-1,1))
diagnosis_cor = data.frame('Module' = gsub('ME','',rownames(moduleTraitCor)),
'MTcor' = moduleTraitCor[,'Diagnosis'],
'MTpval' = moduleTraitPvalue[,'Diagnosis'])
genes_info = genes_info %>% left_join(diagnosis_cor, by='Module')
rm(moduleTraitCor, moduleTraitPvalue, textMatrix, diagnosis_cor, yLabels)Association between Cluster-Diagnosis correlation and percentage of DE genes
**Modules with a high Module-Diagnosis (absolute) correlation should have a high content of differentially expressed genes:
This appears to be the case
plot_data = genes_info %>% group_by(Module, MTcor) %>% summarise(p = 100*mean(significant))
ggplotly(plot_data %>% ggplot(aes(MTcor, p)) + geom_smooth(color='gray', se=FALSE) +
geom_point(color=plot_data$Module, alpha=0.7, aes(id=Module)) +
geom_hline(yintercept=mean(plot_data$p),color='gray',linetype='dotted') +
theme_minimal() +
xlab('Cluster-diagnosis correlation') + ylab('Percentage of DE genes'))3.2.2 Gene-level metrics
Gene significance: is the correlation between the gene and the trait we are interested in. A positive gene significance means the gene is overexpressed and a negative value means its underexpressed. (The term ‘significance’ is not very acurate because it’s not actually measuring statistical significance, it’s just a correlation, but that’s how they call it in WGCNA…)
Module Membership is the correlation of the module’s eigengene and the expression profile of a gene. The higher the Module Membership, the more similar the gene is to the genes that constitute the module.
Note: Some correlations weren’t able to be calculated with the original function (weirdly because the correlation was too strong), so they are calculated again in the end using the polyserial function instead of hetcor, which is robust enough to calculate them
# It's more efficient to iterate the correlations one by one, otherwise it calculates correlations between the eigengenes and also between the genes, which we don't need
# Check if MM information already exists and if not, calculate it
if(file.exists('./../Data/preprocessedData/WGCNA_metrics.csv')){
dataset = read.csv('./../Data/preprocessedData/WGCNA_metrics.csv')
#dataset$Module = dataset[,clustering_selected]
} else {
############# 1. Calculate Gene Significance
GS_info = data.frame('ID'=rownames(datExpr),
'GS'=datExpr %>% apply(1, function(x) hetcor(x,datMeta$Diagnosis)$correlations[1,2])) %>%
mutate('GSpval'=corPvalueStudent(GS, ncol(datExpr)))
############# 2. Calculate Module Membership
#setup parallel backend to use many processors
cores = detectCores()
cl = makeCluster(cores-1)
registerDoParallel(cl)
# Create matrix with MM by gene
MM = foreach(i=1:nrow(datExpr), .combine=rbind) %dopar% {
library(polycor)
tempMatrix = apply(MEs, 2, function(x) hetcor(as.numeric(datExpr[i,]), x)$correlations[1,2])
tempMatrix
}
# Stop clusters
stopCluster(cl)
rownames(MM) = rownames(datExpr)
colnames(MM) = paste0('MM', gsub('ME','',colnames(MEs)))
# Calculate p-values
MMpval = MM %>% corPvalueStudent(ncol(datExpr)) %>% as.data.frame
colnames(MMpval) = paste0('MMpval', gsub('ME','',colnames(MEs)))
MM = MM %>% as.data.frame %>% mutate(ID = rownames(.))
MMpval = MMpval %>% as.data.frame %>% mutate(ID = rownames(.))
# Join and save results
dataset = genes_info %>% dplyr::select(ID, gene.score, Module, module_number, MTcor, MTpval) %>%
left_join(GS_info, by='ID') %>%
left_join(MM, by='ID') %>%
left_join(MMpval, by='ID')
write.csv(dataset, file = './../Data/preprocessedData/WGCNA_metrics.csv', row.names = FALSE)
rm(cores, cl)
}
GS_missing = dataset$ID[is.na(dataset$GS)] %>% as.character
if(length(GS_missing)>0){
cat(paste0(length(GS_missing),' correlations between genes and Diagnosis could not be calculated, ',
'calculating them with the polyserial function'))
for(g in GS_missing){
dataset$GS[dataset$ID == g] = polyserial(as.numeric(datExpr[g,]), datMeta$Diagnosis)
}
}## 2 correlations between genes and Diagnosis could not be calculated, calculating them with the polyserial function## Warning in polyserial(as.numeric(datExpr[g, ]), datMeta$Diagnosis): initial
## correlation inadmissible, -1.00074263127061, set to -0.9999## Warning in polyserial(as.numeric(datExpr[g, ]), datMeta$Diagnosis): initial
## correlation inadmissible, 1.00279583696106, set to 0.9999rm(GS_missing)Association between Gene Significance and cluster-diagnosis correlation
plot_data = dataset %>% dplyr::select(ID, MTcor, GS) %>% mutate(mean_expr = rowMeans(datExpr)) %>%
left_join(genes_info %>% dplyr::select(ID, Group, gene.score), by='ID') %>%
left_join(genes_info %>% dplyr::select(ID, log2FoldChange, shrunken_log2FoldChange, significant,
Module), by='ID') %>%
left_join(data.frame(MTcor=unique(dataset$MTcor)) %>% arrange(by=MTcor) %>%
mutate(order=1:length(unique(dataset$MTcor))), by='MTcor')
plot_data %>% ggplot(aes(MTcor, GS)) + geom_point(color = plot_data$Module, alpha = 0.1) +
geom_smooth(color='gray', se = FALSE) + theme_minimal() +
xlab('Cluster-diagnosis correlation') + ylab('Gene Significance') +
ggtitle(paste0('R^2=',round(cor(plot_data$MTcor, plot_data$GS)^2,2)))
Association to log-fold change
p1 = plot_data %>% ggplot(aes(MTcor, shrunken_log2FoldChange)) + geom_point(color=plot_data$Module, alpha=0.1) +
geom_smooth(color='gray', se=FALSE) + xlab('Cluster-diagnosis correlation') + ylab('Shrunken LFC') +
theme_minimal() + ggtitle(paste0('R^2 = ', round(cor(plot_data$log2FoldChange, plot_data$MTcor)[1]**2,2)))
p2 = plot_data %>% ggplot(aes(GS, shrunken_log2FoldChange)) + geom_point(color=plot_data$Module, alpha=0.1) +
geom_smooth(color='gray', se=FALSE) + xlab('Gene Significance') + ylab('Shrunken LFC') +
theme_minimal() + ggtitle(paste0('R^2 = ', round(cor(plot_data$log2FoldChange, plot_data$GS)[1]**2,2)))
grid.arrange(p1, p2, nrow=1)
rm(p1,p2)Association to mean level of expression
p1 = plot_data %>% ggplot(aes(MTcor, mean_expr)) +
geom_point(color=plot_data$Module, alpha=0.1) + geom_smooth(color='gray', se = TRUE, alpha = 0.2) +
xlab('Cluster-diagnosis correlation') + ylab('Mean level of expression') + theme_minimal() +
ggtitle(paste0('R^2 = ', round(cor(plot_data$log2FoldChange, plot_data$MTcor)[1]**2,2)))
p2 = plot_data %>% ggplot(aes(GS, mean_expr)) +
geom_point(color=plot_data$Module, alpha=0.1) + geom_smooth(color='gray', se = TRUE, alpha = 0.2) +
xlab('Gene Significance') + ylab('Mean level of expression') + theme_minimal() +
ggtitle(paste0('R^2 = ', round(cor(plot_data$log2FoldChange, plot_data$GS)[1]**2,2)))
grid.arrange(p1, p2, nrow=1)
rm(p1,p2)3D plot to try to understand better why the patterns form MTcor and GS with mean expression don’t agree between them (I don’t think this helps)
plot_ly(plot_data, x = ~MTcor, y = ~GS, z = ~mean_expr, color= ~Module, colors = sort(unique(plot_data$Module)),
alpha = 0.5, size = 0.5)## No trace type specified:
## Based on info supplied, a 'scatter3d' trace seems appropriate.
## Read more about this trace type -> https://plot.ly/r/reference/#scatter3d## No scatter3d mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode3.2.3 SFARI Genes
Enrichment in SFARI Genes
The cluster-diagnosis correlation measures the relevance of ASD-related gene expression patterns in a cluster. We now want to use a similar metric to calculate the ``relevance" of SFARI Genes in a cluster.
One way to do this would be to calculate the percentage of genes in each module that belong to the SFARI Genes, but this metric does not take into account the size of the module, so it can favour small clusters, which can have more extreme but less robust results. Instead, the Over Representation Analysis (ORA) does take the size of a cluster into consideration, using the hypergeometric distribution to calculate the probability of a module of size \(n\) containing at least \(s\) SFARI Genes.
To do this, the ORA interprets the number of genes (\(n\)) in a cluster as n random draws without replacement from a finite population of size \(N\), and the number of SFARI genes in the cluster (\(s\)) as s successes in those \(n\) draws, where we know that \(N\) contains exactly \(S\) successes, and it uses the hypergeometric distribution to calculate the statistical significance of having drawn \(s\) successes out of \(n\) draws.
# Calculate % and ORA of SFARI Genes in each module
modules = unique(genes_info$Module[genes_info$Module!='gray']) %>% as.character
# We need the entrez ID of the genes for this
getinfo = c('ensembl_gene_id','entrezgene')
mart = useMart(biomart='ENSEMBL_MART_ENSEMBL', dataset='hsapiens_gene_ensembl',
host='feb2014.archive.ensembl.org')
biomart_output = getBM(attributes=getinfo, filters=c('ensembl_gene_id'),
values=genes_info$ID[genes_info$Module!='gray'], mart=mart) %>%
left_join(genes_info %>% dplyr::select(ID, Module,gene.score), by = c('ensembl_gene_id'='ID'))
# We need to build a term2gene dataframe with the genes and their SFARI Scores
term2gene = biomart_output %>% mutate(term = ifelse(gene.score == 'Others', 'Others', 'SFARI'),
'gene' = entrezgene) %>% dplyr::select(term, gene) %>% distinct
enrichment_data = data.frame('Module' = modules, 'size' = 0, 'perc_SFARI' = 0,
'pval_ORA' = 0, 'padj_ORA' = 0)
for(i in 1:length(modules)){
module = modules[i]
genes_in_module = biomart_output$entrezgene[biomart_output$Module==module]
ORA_module = enricher(gene = genes_in_module, universe = biomart_output$entrezgene %>% as.character,
pAdjustMethod = 'bonferroni', TERM2GENE = term2gene,
pvalueCutoff = 1, qvalueCutoff = 1, maxGSSize = 50000) %>%
data.frame %>% dplyr::select(-geneID,-Description)
ORA_pval = ifelse('SFARI' %in% ORA_module$ID, ORA_module$pvalue[ORA_module$ID=='SFARI'], 1)
ORA_padj = ifelse('SFARI' %in% ORA_module$ID, ORA_module$p.adjust[ORA_module$ID=='SFARI'], 1)
enrichment_data[i,-1] = c(length(genes_in_module),
mean(genes_info$gene.score[genes_info$Module==module]!='Others'),
ORA_pval, ORA_padj)
}
enrichment_data = enrichment_data %>%
left_join(genes_info %>% dplyr::select(Module, MTcor) %>% unique, by = 'Module')
write.csv(enrichment_data, file = './../Data/preprocessedData/SFARI_enrichment_by_cluster.csv')
rm(i, module, genes_in_module, ORA_module, ORA_pval, ORA_padj, getinfo, mart, term2gene)ggplotly(enrichment_data %>% ggplot(aes(perc_SFARI, 1-pval_ORA, size=size)) +
geom_point(color = enrichment_data$Module, alpha = 0.5, aes(id=Module)) +
geom_smooth(color='gray', se = TRUE, alpha = 0.1) + coord_cartesian(ylim=c(0,1)) +
xlab('Percentage of SFARI Genes') + ylab('ORA Enrichment') + theme_minimal() +
theme(legend.position = 'none') + ggtitle('Modules'))## Warning: Ignoring unknown aesthetics: id## `geom_smooth()` using method = 'loess' and formula 'y ~ x'Association between the Enrichment in SFARI Genes and cluster-diagnosis correlation
There doesn’t seem to be a strong relation between these two metrics, if anything, they seem to contradict each other, mainly in modules with a high Module-Diagnosis correlation, which are the ones with the lowest trend for Enrichment on SFARI Genes.
ggplotly(enrichment_data %>% ggplot(aes(MTcor, 1-pval_ORA, size=size)) +
geom_point(color=enrichment_data$Module, alpha=0.5, aes(id=Module)) +
geom_smooth(color='#cccccc', size = 0.5,alpha=0.1) +
xlab('Module-Diagnosis Correlation') + ylab('Enrichment in SFARI Genes') +
theme_minimal() + theme(legend.position = 'none') + coord_cartesian(ylim=c(0,1)) +
ggtitle('SFARI Enrichment vs Module-Diagnosis Correlation'))Association between the Enrichment in SFARI Genes and the mean expression of the cluster
There is a positive relation between these two metrics; the higher the mean expression of the genes in a module, the higher the enrichment in SFARI Genes, especially for modules with medium to high level of expression.
mean_expr_by_module = data.frame(ID = rownames(datExpr), mean_expr = rowMeans(datExpr)) %>%
left_join(dataset %>% dplyr::select(ID, Module), by = 'ID') %>%
group_by(Module) %>% summarise(mean_expr = mean(mean_expr))
ggplotly(enrichment_data %>% left_join(mean_expr_by_module, by = 'Module') %>%
ggplot(aes(mean_expr, 1-pval_ORA, size=size)) +
geom_point(color=enrichment_data$Module, alpha=0.5, aes(id=Module)) +
geom_smooth(color='#cccccc', size = 0.5,alpha=0.1) + coord_cartesian(ylim=c(0,1)) +
xlab('Mean level of expression') + ylab('Enrichment in SFARI Genes') +
theme_minimal() + theme(legend.position = 'none') +
ggtitle('SFARI Enrichment vs Module-Diagnosis Correlation'))Association between WGCNA metrics and SFARI Genes
wt = plot_data %>% mutate(abs_GS = abs(GS), Group = factor(Group, levels = c('SFARI','Neuronal','Others'))) %>%
wilcox_test(abs_GS~Group, p.adjust.method='BH') %>% add_x_position(x = 'group')
increase = 0.1
base = 1.05
pos_y_comparisons = c(base, base + increase, base)
p1 = plot_data %>% mutate(Group = factor(Group, levels = c('SFARI', 'Neuronal', 'Others'))) %>%
ggplot(aes(Group, abs(GS))) +
geom_boxplot(outlier.colour='#cccccc', outlier.shape='o', outlier.size=3, aes(fill=Group)) +
stat_pvalue_manual(wt, y.position = pos_y_comparisons, tip.length = .02) +
scale_fill_manual(values=c('#00A4F7', SFARI_colour_hue(r=c(8,7)))) + theme_minimal() +
scale_x_discrete(labels=c('SFARI' = 'SFARI\nGenes', 'Others' = 'Other\ngenes',
'Neuronal' = 'Neuronal\ngenes')) +
ylab('Gene Significance magnitude') + xlab('SFARI Scores') + theme(legend.position='none')
wt = plot_data %>% mutate(abs_MTcor=abs(MTcor), Group=factor(Group, levels=c('SFARI','Neuronal','Others'))) %>%
wilcox_test(abs_MTcor~Group, p.adjust.method='BH') %>% add_x_position(x = 'group')
increase = 0.08
base = 1.05
pos_y_comparisons = c(base, base + increase, base)
p2 = plot_data %>% mutate(Group = factor(Group, levels = c('SFARI', 'Neuronal', 'Others'))) %>%
ggplot(aes(Group, abs(MTcor))) +
geom_boxplot(outlier.colour='#cccccc', outlier.shape='o', outlier.size=3, aes(fill=Group)) +
stat_pvalue_manual(wt, y.position = pos_y_comparisons, tip.length = .02) +
scale_fill_manual(values=c('#00A4F7', SFARI_colour_hue(r=c(8,7)))) + theme_minimal() +
scale_x_discrete(labels=c('SFARI' = 'SFARI\nGenes', 'Others' = 'Other\ngenes',
'Neuronal' = 'Neuronal\ngenes')) +
ylab('Cluster-diagnosis correlation magnitude') + xlab('SFARI Scores') + theme(legend.position='none')
grid.arrange(p2, p1, nrow=1)
rm(p1, p2, increase, base, pos_y_comparisons, wt)Association between WGCNA metrics and SFARI Scores
comparisons = list(c('1','2'), c('2','3'), c('3','Neuronal'), c('Neuronal','Others'),
c('1','3'), c('3','Others'), c('2','Neuronal'),
c('1','Neuronal'), c('2','Others'), c('1','Others'))
increase = 0.1
base = 1.05
pos_y_comparisons = c(rep(base, 4), rep(base + increase, 2), base + 2:5*increase)
p1 = plot_data %>% ggplot(aes(gene.score, abs(GS), fill=gene.score)) +
geom_boxplot(outlier.colour='#cccccc', outlier.shape='o', outlier.size=3) +
stat_compare_means(comparisons = comparisons, label = 'p.signif', method = 't.test',
method.args = list(var.equal = FALSE), label.y = pos_y_comparisons,
tip.length = .02) +
scale_fill_manual(values=SFARI_colour_hue(r=c(1:3,8,7))) + theme_minimal() +
scale_x_discrete(labels=c('1' = 'SFARI\nScore 1', '2' = 'SFARI\nScore 2', '3' = 'SFARI\nScore 3',
'Others' = 'Other\ngenes', 'Neuronal' = 'Neuronal\ngenes')) +
ylab('Gene Significance magnitude') + xlab('SFARI Scores') + theme(legend.position='none')
increase = 0.08
base = 1
pos_y_comparisons = c(rep(base, 4), rep(base + increase, 2), base + 2:5*increase)
p2 = plot_data %>% ggplot(aes(gene.score, abs(MTcor), fill=gene.score)) +
geom_boxplot(outlier.colour='#cccccc', outlier.shape='o', outlier.size=3) +
stat_compare_means(comparisons = comparisons, label = 'p.signif', method = 't.test',
method.args = list(var.equal = FALSE), label.y = pos_y_comparisons,
tip.length = .01) +
scale_fill_manual(values=SFARI_colour_hue(r=c(1:3,8,7))) + theme_minimal() +
scale_x_discrete(labels=c('1' = 'SFARI\nScore 1', '2' = 'SFARI\nScore 2', '3' = 'SFARI\nScore 3',
'Others' = 'Other\ngenes', 'Neuronal' = 'Neuronal\ngenes')) +
ylab('Cluster-diagnosis correlation magnitude') + xlab('SFARI Scores') + theme(legend.position='none')
grid.arrange(p2, p1, nrow=1)
rm(p1, p2, comparisons, increase, base, pos_y_comparisons)Session info
sessionInfo()## R version 3.6.3 (2020-02-29)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.6 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
##
## locale:
## [1] LC_CTYPE=en_GB.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB.UTF-8 LC_COLLATE=en_GB.UTF-8
## [5] LC_MONETARY=en_GB.UTF-8 LC_MESSAGES=en_GB.UTF-8
## [7] LC_PAPER=en_GB.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] doParallel_1.0.15 iterators_1.0.13 foreach_1.5.1
## [4] clusterProfiler_3.12.0 biomaRt_2.40.5 rstatix_0.6.0
## [7] polycor_0.7-10 expss_0.10.7 WGCNA_1.69
## [10] fastcluster_1.2.3 dynamicTreeCut_1.63-1 ggExtra_0.9
## [13] ggpubr_0.2.5 magrittr_2.0.1 GGally_1.5.0
## [16] colorspace_2.0-2 gridExtra_2.3 viridis_0.6.1
## [19] viridisLite_0.4.0 RColorBrewer_1.1-2 dendextend_1.15.1
## [22] plotly_4.9.2 glue_1.4.2 reshape2_1.4.4
## [25] forcats_0.5.0 stringr_1.4.0 dplyr_1.0.1
## [28] purrr_0.3.4 readr_1.3.1 tidyr_1.1.0
## [31] tibble_3.1.2 ggplot2_3.3.5 tidyverse_1.3.0
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.2 tidyselect_1.1.1
## [3] RSQLite_2.2.0 AnnotationDbi_1.46.1
## [5] htmlwidgets_1.5.3 grid_3.6.3
## [7] BiocParallel_1.18.1 munsell_0.5.0
## [9] codetools_0.2-16 preprocessCore_1.46.0
## [11] miniUI_0.1.1.1 withr_2.4.2
## [13] GOSemSim_2.10.0 Biobase_2.44.0
## [15] highr_0.9 knitr_1.32
## [17] rstudioapi_0.13 stats4_3.6.3
## [19] ggsignif_0.6.2 DOSE_3.10.2
## [21] labeling_0.4.2 urltools_1.7.3
## [23] GenomeInfoDbData_1.2.1 polyclip_1.10-0
## [25] bit64_4.0.5 farver_2.1.0
## [27] vctrs_0.3.8 generics_0.1.0
## [29] xfun_0.25 GenomeInfoDb_1.20.0
## [31] R6_2.5.1 graphlayouts_0.7.0
## [33] locfit_1.5-9.4 DelayedArray_0.10.0
## [35] bitops_1.0-7 cachem_1.0.6
## [37] reshape_0.8.8 fgsea_1.10.1
## [39] gridGraphics_0.5-1 assertthat_0.2.1
## [41] promises_1.2.0.1 scales_1.1.1
## [43] ggraph_2.0.3 nnet_7.3-14
## [45] enrichplot_1.4.0 gtable_0.3.0
## [47] tidygraph_1.2.0 rlang_0.4.11
## [49] genefilter_1.66.0 splines_3.6.3
## [51] lazyeval_0.2.2 acepack_1.4.1
## [53] impute_1.58.0 broom_0.7.0
## [55] europepmc_0.4 checkmate_2.0.0
## [57] BiocManager_1.30.16 yaml_2.2.1
## [59] abind_1.4-5 modelr_0.1.6
## [61] crosstalk_1.1.1 backports_1.2.1
## [63] httpuv_1.6.1 qvalue_2.16.0
## [65] Hmisc_4.4-0 tools_3.6.3
## [67] ggplotify_0.1.0 ellipsis_0.3.2
## [69] jquerylib_0.1.4 BiocGenerics_0.30.0
## [71] ggridges_0.5.3 Rcpp_1.0.7
## [73] plyr_1.8.6 zlibbioc_1.30.0
## [75] base64enc_0.1-3 progress_1.2.2
## [77] RCurl_1.98-1.4 prettyunits_1.1.1
## [79] rpart_4.1-15 cowplot_1.1.1
## [81] S4Vectors_0.22.1 SummarizedExperiment_1.14.1
## [83] haven_2.2.0 ggrepel_0.9.1
## [85] cluster_2.1.0 fs_1.5.0
## [87] data.table_1.14.0 DO.db_2.9
## [89] openxlsx_4.2.4 triebeard_0.3.0
## [91] reprex_0.3.0 mvtnorm_1.1-2
## [93] matrixStats_0.60.1 hms_1.1.0
## [95] mime_0.11 evaluate_0.14
## [97] xtable_1.8-4 XML_3.99-0.3
## [99] rio_0.5.16 jpeg_0.1-9
## [101] readxl_1.3.1 IRanges_2.18.3
## [103] compiler_3.6.3 crayon_1.4.1
## [105] htmltools_0.5.2 mgcv_1.8-36
## [107] later_1.3.0 Formula_1.2-4
## [109] geneplotter_1.62.0 lubridate_1.7.10
## [111] DBI_1.1.1 tweenr_1.0.2
## [113] dbplyr_1.4.2 MASS_7.3-53
## [115] Matrix_1.2-18 car_3.0-7
## [117] cli_3.0.1 igraph_1.2.6
## [119] GenomicRanges_1.36.1 pkgconfig_2.0.3
## [121] rvcheck_0.1.8 foreign_0.8-76
## [123] xml2_1.3.2 annotate_1.62.0
## [125] bslib_0.3.0 XVector_0.24.0
## [127] rvest_0.3.5 yulab.utils_0.0.2
## [129] digest_0.6.27 rmarkdown_2.7
## [131] cellranger_1.1.0 fastmatch_1.1-3
## [133] htmlTable_1.13.3 curl_4.3.2
## [135] shiny_1.6.0 nlme_3.1-153
## [137] lifecycle_1.0.0 jsonlite_1.7.2
## [139] carData_3.0-4 fansi_0.5.0
## [141] pillar_1.6.2 lattice_0.20-41
## [143] fastmap_1.1.0 httr_1.4.2
## [145] survival_3.2-7 GO.db_3.8.2
## [147] zip_2.2.0 UpSetR_1.4.0
## [149] png_0.1-7 bit_4.0.4
## [151] ggforce_0.3.1 stringi_1.7.4
## [153] sass_0.4.0 blob_1.2.2
## [155] DESeq2_1.24.0 latticeExtra_0.6-29
## [157] memoise_2.0.0